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Abstract. We calculate the redshift-space power spectrum of the Sloan Digital Sky Survey (SDSS) Data Release 4 (DR4) 
Luminous Red Galaxy (LRG) sample, finding evidence for a full series of acoustic features down to the scales of ~ 0.2 h Mpc“*. 
This corresponds up to the 7th peak in the CMB angular power spectrum. The acoustic scale derived, (105.4 ± 2.3) Mpc, 
agrees very well with the “concordance” model prediction and also with the one determined via the analysis of the spatial 
two-point correlation function by Eisenstein et al. (2005). The models with baryonic features are favored by 3.3(T over their 
“smoothed-out” counterparts without any oscillatory behavior. This is not only an independent confirmation of Eisenstein et al. 
(2005) results made with different methods and software but also, according to our knowledge, the first determination of the 
power spectrum of the SDSS LRG sample. 
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1. Introduction 

In the beginning of 1970’s it was already realized that acoustic 
waves in the tightly coupled baryon-photon fluid prior to the 
epoch of recombination will lead to the characteristic maxima 
and minima in the post-recombination matter power spectrum. 
The same mechanism is also responsible for the prominent 
peak structure in the CMB angular power spectrum (Sunyaev 
& Zeldovich 1970; Peebles & Yu 1970; Doroshkevich et al. 
1978). The scale of these features reflects the size of the sound 
horizon, which itself is fully determined given the physical den¬ 
sities Q.bh^ and The acoustic horizon can be calibrated 

using the CMB data, thus turning it into a standard ruler which 
can be used to carry out classical cosmological tests. For exam¬ 
ple, if we are able to measure the redshift and angular intervals 
corresponding to the physically known acoustic scale in the 
matter power spectrum at a range of redshifts, we can immedi¬ 
ately And angular diameter distance dp, and Hubble parameter 
// as a function of redshift. Having good knowledge of these 
dependencies allows us to put constraints on the properties of 
the dark energy. To carry out this project one needs a tracer 
population of objects whose clustering properties with respect 
to the underlying matter distribution is reasonably well un¬ 
derstood. There have been several works discussing the usage 
of galaxies (Blake & Glazebrook 2003; Hu & Haiman 2003; 
Linder 2003; Seo & Eisenstein 2003) and clusters of galaxies 
(Hu & Haiman 2003; Majumdar & Mohr 2004; Hiitsi 2005) 
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for this purpose. What is most important is that already cur¬ 
rently existing galaxy redshift surveys have lead to the detec¬ 
tion of acoustic features in the spatial distribution of galaxies, 
this way providing clearest support for the standard gravita¬ 
tional instability picture of the cosmic structure formation. In 
the paper by Eisenstein et al. (2005) the detection of the acous¬ 
tic “bump” in the two-point redshift-space correlation function 
of the SDSS ' LRG sample is announced. The discovery of 
similar features in the power spectrum of 2dE ^ galaxies is 
presented in Cole et al. (2005). These results clearly demon¬ 
strate the great promise of the future dedicated galaxy redshift 
surveys like K.A.O.S.^ Similarly, useful measurements of the 
acoustic scale can be hoped by the planned SZ cluster surveys 
like the ones carried out by the PLANCK Surveyor space¬ 
craft and SPT ^ (Hiitsi 2005) and also with a large future pho¬ 
tometric redshift surveys (Blake & Bridle 2005). Eor the SZ 
surveys one needs an additional optical follow-up to get esti¬ 
mates for the cluster redshifts. In this paper we calculate the 
redshift-space power spectrum of the SDSS LRG sample find¬ 
ing evidence for the acoustic oscillations down to the scales of 
~ 0.2/iMpc *, which effectively correspond up to the 7. peak 
in the CMB angular power spectrum. These scales in the CMB 
are very strongly damped due to the finite width of the last¬ 
scattering surface and also due to the Silk damping (Silk 1968). 

' http://www.sdss.org/ 

^ http://www.mso.anu.edu.au/2dFGRS/ 

^ http://www.noao.edu/kaos/ 
http://astro.estec.esa.nl/Planck 
^ http://astro.uchicago.edu/spt 
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Fig. 1 Acoustic oscillations in the CMB (upper panel) and linear mat¬ 
ter power spectrum (lower panel) for the “concordance” cosmological 
model. Here, as we have plotted the spectra against spatial wavenum¬ 
ber k, we have changed the standard notation of C{ to Cj. Due to the 
factor the first CMB acoustic peak is barely visible. Density fluctua¬ 
tions in matter at smaller scales, being mostly induced by the velocity 
fields, are out of phase with respect to the fluctuations in the CMB 
component. Also the fluctuation period is twice as large. 


This can be seen in Fig. 1 ® where the CMB data is plotted in a 
somewhat unusual way to enhance the acoustic features at the 
high wavenumber damping tail. Also, at those scales the sec¬ 
ondary CMB anisotropies (mostly thermal Sunyaev-Zeldovich 
effect (Sunyaev & Zeldovich 1972, 1980)) start to dominate 
over the primary signal. On the other hand, features in the mat¬ 
ter power spectrum, although being small (~ 5% fluctuations), 
are preserved by the linear evolution and so opening up the way 
to probe acoustic phenomena at scales smaller than the ones ac¬ 
cessible for the CMB studies. 

The paper is structured as follows. In Sec. 2 we describe the 
dataset to be analyzed. Sec. 3 presents the method of the power 
spectrum calculation. In Sec. 4 we determine power spectrum 
errors and covariance matrix. Sec. 5 discusses the convolution 
effect of the survey window. Analytical model spectra are pre¬ 
sented in Sec. 6. The results of the measurement of the acoustic 
scale are given in Sec. 7. Correlation function analysis is car¬ 
ried out in Sec. 8. In Sec. 9 we compare the measured power 
spectrum with the published results for the 2dF and SDSS main 
sample, and finally we conclude with Sec. 10. 

2. Data 

We analyze the publicly available data from the SDSS DR4 
(Adelman-McCarthy et al. 2005). Specifically, we carry out our 
power spectrum measurements using the subset of the SDSS 
spectroscopic sample known as the Luminous Red Galaxy 
(LRG) sample. The LRG selection algorithm (Eisenstein et al. 
2001) selects ~ 12 galaxies per square degree meeting spe- 

® Here instead of the usual multipole number ( we have plotted 
the CMB angular power spectrum against the wavenumber k. For the 
“concordance” cosmological model ( = 9990A:[/!Mpc“*]. 



Comoving Distance [h'^Mpc] 

Fig. 2 Comoving number density of galaxies as a function of comov¬ 
ing distance. Smooth solid line shows a cubic spline fit to the number 
density estimated for 50 discrete radial bins. 

cific color and magnitude criteria The resulting set of galax¬ 
ies consists mostly of an early types populating dense cluster 
environments and as such are significantly biased (bias factor 
b ~ 2) with respect to the underlying matter distribution. The 
selection method is very effective producing a galaxy sample 
with a reasonably high density up to the redshift of z ~ 0.5. 

Since the selection criteria are very complicated, involving 
both cuts in magnitude and in color, and also due to the steep¬ 
ness of the luminosity function the usual method of using only 
the luminosity function to determine radial selection function 
does not work here (Zehavi et al. 2005). Here we simply build 
the radial selection function as a smooth spline fit to the num¬ 
ber density profiles shown in Fig. 2. To calculate distances we 
choose the cosmological parameters as given by the WMAP * 
“concordance” model (Spergel et al. 2003). Unfortunately the 
coverage masks of the SDSS DR4 spectroscopic sample are 
not available in a readily accessible format and so we chose to 
build the angular survey masks using the galaxy data itself 
As the number density of galaxies in the sample is rather high, 
one can determine relatively accurately the beginning, ending 
and also possible gaps in the scan stripes. We have built angu¬ 
lar masks using both the whole DR4 galaxy sample and LRGs 
only. The measured power spectra are practically identical with 
only some minor differences on smaller scales (see Fig. 6). This 
can be seen as an indication that our power spectrum measure¬ 
ments are rather stable against small uncertainties in the survey 
geometry. The angular distribution of the galaxies and also the 

^ For the exact details of the selection criteria see Eisenstein et al. 
( 2001 ) 

* http://lambda.gsfc.nasa.gov/product/map/ 

* In principle one can build the angular masks using the raw tiling 
information, but as we show later our approximate treatment is prob¬ 
ably rather fine, since the results seem to be quite stable against small 
uncertainties in the mask. More rigorous approach should certainly 
address the issues of survey boundaries and completeness fluctuations 
(expected to be small due to the very effective tiling algorithm by 
Blanton et al. (2003)) in a much better detail. 
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Fig. 3 Angular distribution of galaxies given in the SDSS survey coor¬ 
dinates (d, 7/). The survey mask is shown with solid lines. The vertical 
dashed lines show the division of the sample into 22 separate regions 
each containing ~ 2350 galaxies. This division will be exploited in the 
“jackknife” error analysis of the correlation function. 


boundaries of the survey mask built in the above mentioned 
way (here using all the galaxies) is shown in Fig. 3. Here the 
angular positions are plotted using the so-called survey coordi¬ 
nate system of the SDSS 

We have selected all the objects that have spectrum classi- 
hed as galaxy (i.e. SpecClass=2) and are additionally flagged 
as GALAXY_RED or GALAXY_RED_II (i.e. PrimTarget bit mask 
set as ®x2® or ®x4®®®®®®, respectively). Only galaxies for 
which the redshift conhdence parameter, zConf, is greater than 
0.95 were used. We apply lower and upper redshift cutoffs of 
0.16 and 0.47 as also done in Eisenstein et al. (2005). The lower 
cutoff is needed since the color cuts that dehne the LRG sample 
break down for redshifts below ~ 0.2 (Eisenstein et al. 2001). 
Eor the analysis presented in this paper we have excluded the 
three southern stripes since these just increase the sidelobes of 
the survey window without adding much of the extra volume. 
We have also removed some minor parts of the sample to ob¬ 
tain more continuous and smooth chunk of volume. In total the 
analyzed galaxy sample covers ~ 0.75 hT^ Gpc^ over ~ 3850 
square degrees on the sky and contains 51,763 galaxies. 


3. Power spectrum calculation 

We calculate the power spectrum using a direct Eourier method 
as described in Eeldman et al. (1994) (EKP). Strictly speak¬ 
ing, power spectra determined this way are the so-called 
pseudospectra, meaning that the estimates derived are con¬ 
volved with a survey window. Since in the case of the ana¬ 
lyzed LRG sample the volume covered is very large, reaching 
0.75 Gpc^, and also the survey volume has relatively large 
dimensions along all perpendicular directions, the correlations 
in the Eourier space are rather compact. On intermediate scales 
and in the case the power spectrum binning is chosen wide 

The transformations between various coordinate systems used by 
the SDSS are given e.g. in Stoughton et al. (2002). 


enough, EKP estimator gives a good approximation to the true 
underlying power. 

The EKP estimate for a 3D pseudospectrum reads as: 


P(k) = |P(k)|2-P,hot, 

(1) 

where 


F{k) - J" d^rF(r) exp(/k ■ r). 

(2) 

Here F(r) is the weighted density contrast field: 


F(r) = w(r) [ng(r) - a«,(r)] . 

(3) 

ng(r) and ni(r) denote the number densities of the 
galaxy catalog and a synthetic random catalog with 
selection criteria, respectively. Since we are dealing 
Crete point processes, densities can be given as: 

analyzed 
the same 

with dis- 

= X ’ 

i 

(4) 

n,(r) = 2Ar-<), 

(5) 


where and r* denote the location of the /-th point in real and 
synthetic catalog, respectively, and 6^ is the 3D Dirac delta 
function, a in Eq. (3) is the ratio of the number of galaxies to 
the number of random points in the synthetic catalog i.e. a - 
In our calculations we have Ns - 10 and thus a ^ 0.0052. 
Eor the weight function w(r) there have been traditionally three 
choices in the literature: 

for volume weighting 
for number weighting (6) 

for an optimal EKP weighting. 

Here h(r) is the average number density of galaxies at comov¬ 
ing location r i.e. the radial selection function of the survey (see 
Pig. 2) times the angular mask (Pig. 3). In our calculations we 
use an optimal EKP weighting scheme, although pure volume 
weighting would give practically the same results, especially 
on the larger scales {k < 0.09 h Mpc ’), since then for the ma¬ 
jority of the sample n(r)P ~ 3 ". The weights in Eq. (3) are 
normalized such that: 

dVh2(r)wV) = 1, (7) 

which can be approximated as the following sum over the syn¬ 
thetic catalog 

^ I . ( 8 ) 

i 

The last term in Eq. (1) represents the Poissonian discreteness 
noise and can be expressed as: 

Pshot = (1 + O') J" dV«(r)w^(r) a(l -i- a) ^ w'^(rj). (9) 

Including all the modes down to the scales of fc ~ 0.25 AMpc^' 
the effective value for h(r)P drops down to ~ 1.5. 

We assume that the survey selection does not have any other angu¬ 
lar dependence except for the applied angular mask i.e. we can replace 
r) by the modulus rj. 


h(r) 

w(r) oc •; const 


l+n(r)/> 









4 


G. Hiitsi: SDSS LRG acoustic oscillations 


Since we are using Fast Fourier Transforms (FFTs) to speed 
up the calculation of the Fourier sums, we have to deal with 
some extra complications. As the density field is now “re¬ 
stricted to live” on a regular grid with a finite cell size, we 
have to correct for the smoothing effect this has caused. Also, if 
our underlying density field contains spatial modes with higher 
frequency than our grid’s Nyquist frequency, ki^y, then these 
will be “folded back” into the frequency interval the grid can 
support, increasing power close to ^Ny- the so-called aliasing 
effect. The relation between the spectra calculated using di¬ 
rect summation and the ones found using FFT techniques was 
worked out by Jing (2005). It can be expressed as follows: 

l^(k)lFFT = X + 2kNyn)|2p(k + 2kNyn) + 

neZ 

f^shotXl'^^k + 2kNyn)|2, (10) 

neZ 

where ^(k) is the mass assignment function used to build den¬ 
sity grid out of the point set. We use the Triangular Shaped 
Cloud (TSC) assignment method (Hockney & Eastwood 1988). 
Since the TSC filter can be obtained by convolving uniform 
cube (the Nearest Grid Point filter) two times with itself, the 
Fourier representation of it follows immediately: 


1. Determination of the survey selection function i.e. mean 
underlying number density h(r) (including the survey ge¬ 
ometry), 

2. Calculation of the overdensity field on a grid using TSC 
mass assignment scheme, 

3. Fourier transformation of the gridded density field, 

4. Calculation of the raw 3D power spectrum |F'(k)|ppj,, 

5. Subtraction of the shot noise component from the raw spec¬ 
trum, 

6. Recovery of the angle averaged pseudospectrum P{k) using 
an iterative method of Jing (2005). 

We have applied the above described power spectrum cal¬ 
culation method to a multitude of test problems, the results 
of which can be found in Hiitsi (2005). In Appendix A we 
show only one example, where we successfully recover the un¬ 
derlying power spectrum of galaxy clusters from the VIRGO 
Hubble Volume simulations after applying the selection cri¬ 
teria given in Figs. 2 and 3. 

4. Power spectrum errors and covariance matrix 

We determine power spectrum errors by three different meth¬ 
ods: 


"W(k) = 



k = {ky,k2,h). 


1. Prescription given by FKP that assumes the underlying den¬ 
sity field to be Gaussian. This method also does not treat 
redshift space distortions. Under those simplifying assump¬ 
tions the power spectrum variance can be expressed as: 


Here the sum that represents the contribution from aliases runs 
over all the integer vectors n. Eq. (10) is the direct analog of 
the previous Eq. (1). The convolution with the mass assignment 
filter has introduced factors both to the power spectrum 

and to the shot noise term. The sum in the last term of Eq. (10) 
can be performed analytically for the TSC filter to yield the 
result (Jing 2005) 



= — Xi Z - k") + 5(k' - k")P, (13) 


k k' k" 




2(k) = I d^rn^(r)w^(r)exp(/k ■ r) ^ 


a Z n{r])v?{r]) exp(/k ■ rj), (14) 


5(k) = (1 -I- a) I d^r n(r)w^(r) exp(/k ■ r) ^ 


a(l -H a) ^ exp(/k ■ rj). (15) 


To recover the angle averaged pseudospectrum P{k) from Eq. 
(10) we use an iterative scheme as described in Jing (2005) with 
a slight modification: we do not approximate the small scale 
spectrum by a simple power law, but also allow for the possible 
running of the spectral index i.e. the parametric shape of the 
power spectrum is taken to be a parabola in log-log. Since on 
small scales the power spectrum is dropping fast, the sum over 
n in Eq. (10) is converging rather rapidly. In calculations we 
use only integer vectors with |n| < 5. The angular average is 
taken over all the vectors k laying in the same k-space shell 
with width Ak. The resulting P is taken to be an estimate for 
the pseudospectrum at the wavenumber keff that corresponds to 
the average length of the k-vectors in that shell. 

To summarize, our power spectrum calculation consists of 
the following steps: 

For the NGP filter this sum equals 1, and so one recovers the orig¬ 
inal shot noise term in Eq.(l). 


Here the sum is over all the wavevectors k' and k" pop¬ 
ulating the same k-space shell with radius k and thick¬ 
ness Ak, and Nk denotes the total number of modes in 
that shell. Since the direct summation over all the vec¬ 
tor pairs k' and k" is very slow for the wide k-space 
shells and 512^ grid we use, a Monte Carlo sum is per¬ 
formed instead. Thus we calculate the average of the quan¬ 
tity |P(k)2(k' - k") -H 5 (k' - k")P over the random pairs 
of vectors k' and k" from the same shell. Eor the result to 
converge properly we need on average ~ 10^ random pairs. 

2. The second method is a simple analytical approximation to 
the first one, also due to EKP (see also Tegmark et al. 1998). 
Here the variance is given as: 



2P\k) 

VeffVk ’ 


(16) 


http://www.mpa-garching.mpg.de/Virgo/ 
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Fig. 4 Power spectrum of the SDSS LRG sample with the bin width 
AA: =: 0.005 h Mpc^'. The upper solid line shows the best fitting model 
spectrum and the lower one corresponds to the linearly evolved matter 
power spectrum of the “concordance” cosmological model multiplied 
by the square of the bias parameter b = 1.95. Both of the spectra 
are convolved with a survey window. The dashed lines represent the 
corresponding unconvolved spectra. 



Fig. 5 Power spectrum of the SDSS LRG sample with the bin width 
Afc =: 0.02/iMpc^'. The upper solid line shows the best fitting model 
spectrum and the lower one corresponds to the linearly evolved matter 
power spectrum of the “concordance” cosmological model multiplied 
by the square of the bias parameter b = 1.95. Both of the spectra 
are convolved with a survey window. The dashed lines represent the 
“smoothed-out” versions of the above model spectra. The dotted line 
is the cubic spline fit to the data points. 


where - Ank^Ak/{2 ji)^ is the volume of the At-space shell 
and Veff is the effective volume given by; 

[ f d^rn^(r)w^(r)l 

Veff = -^^- ^2 . (17) 

/dVn4(r)w4(r)[l + ^] 

3. The third method is a Monte Carlo approach that uses 1000 
mock catalogs generated in the way described in Appendix 
B. Here, as we use the 2nd order Lagrangian perturba¬ 
tion theory, we get a good approximation for the mode¬ 
mode couplings that are induced during the quasi-nonlinear 
regime of the evolution of the density fluctuations. Also 
the large-scale redshift distortions are properly accounted 
for. In terms of the Halo Model (see Appendix C) we can 
say that halo-halo clustering term is relatively well approx¬ 
imated. Contributions from the one-halo term can be added 
later, as these allow an analytic treatment. 

The results of the power spectra for the SDSS DR4 LRG 
sample are shown in Figs. 4 and 5. In Fig. 4 the bin width 
Ak ^ 0.005/!Mpc-\ while in Fig. 5 Ak ^ 0.02/iMpc^'. With 
different lines we have shown various model spectra, which 
will be the topic of Sec. 6. 

The comparison of the power spectrum errorbars calculated 
in the above described different ways is provided in Fig. 6. We 
see that the various error estimates are in a very good agree¬ 
ment. In the following we will use only the errorbars given by 
the 3rd method. 

So far we have only found the diagonal terms of the covari¬ 
ance matrix. In order to answer the question of how strongly 
different power spectrum bins are correlated, we have to go a 
step further, and try to estimate the full covariance matrix. 



Fig. 6 The comparison of the different power spectrum error estimates. 
For clarity slight relative shifts of the data points have been applied. 
The errorbars resulting from the 1st method are the rightmost ones and 
the ones from the 3rd method are displayed in the middle. The lines 
show a cubic spline fits to the data points. The solid line corresponds to 
the case when all the available galaxy data is used to find the angular 
mask of the survey, while the dashed line represents the case when 
LRGs only are used for this purpose. 


The FKP result for the full covariance matrix, C,y, 
pie generalization of the Eq. (13): 


Co = 




■zz 




Iki + k. 


2(k'-k") + 5(k'-k") 


,(18) 


where the k-vectors k' and k" lie in shells with width Ak and 
radii k, and kj, respectively. The FKP approach, as mentioned 
above, does not treat mode couplings arising from the nonlinear 






















































6 


G. Hiitsi: SDSS LRG acoustic oscillations 




10 “ 


10 ^ 

0.9 

10 “ 

0.8 


0.7 

10 

0.6 

10 ^ 

0.5 

io “ 

0.4 


0.3 


0.2 


0.1 



1 

0.9 

0.8 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 

0 



0.99995 

0.9999 

0.99985 

0.9998 

0.99975 

0.9997 

0.99965 

0.9996 

0.99955 



0.99995 

0.9999 

0.99985 

0.9998 

0.99975 

0.9997 

0.99965 

0.9996 

0.99955 


Fig. 7 Covariance (left column) and correlation matrices (right column). Top row represents the results from FKP prescription (see Eq. (18)) 
and the middle row the ones from 1000 mock catalogs. The last row displays the nonlinear contribution due to the 1-halo term. 


evolution and also from the redshift space distortions. Linear 
redshift space distortions can be, in principle, included into the 
FKP estimate for the covariance matrix. One can generalize 
the results presented in the Appendix of Zaroubi & Hoffman 
(1996), where the covariance matrix for the Fourier modes 
has been found. Since linear redshift distortions applied on a 
Gaussian field do not change the Gaussianity property, one can 
still use the result from the Appendix B of FKP that relates the 


power spectrum covariance matrix to the covariance matrix of 
the Fourier modes. Also one has to add the shot noise terms 
to the result of Zaroubi & Hoffman (1996). We have carried 
out this exercise, leading us to the high dimensional integrals 
(up to 12 dim.) that turn out to be too time consuming to solve 
in practice. As from the mock catalogs we can hopefully obtain 
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more realistic estimate for the covariance matrix we have not 
followed this path any further. 

The results for the covariance matrix calculation are given 
in Fig. 7. Here the left hand column shows the covariance and 
the right hand column respective correlation matrices: 


Ci. 


y[Cu^j 


(19) 


The power spectrum binning is the same as shown in Fig. 5 i.e. 
Ak ^ 0.02/iMpc *. The top row represents the results from 
Eq. (18), while the middle row the ones from mock catalogs. 
Although the diagonal terms of the covariance matrices in the 
1st and 2nd row are in a very good agreement (see Fig. 6), the 
off-diagonal components differ strongly. This can be explained 
as the result of the extra mode-mode couplings that are not ac¬ 
counted for by the FKP approach. We see that even well sepa¬ 
rated power spectrum bins can be correlated at 30 ... 40% level. 
The bottommost row in Fig. 7 represents the nonlinear contri¬ 
bution to the covariance matrix arising from the 1-halo term 
(see Appendix C). We see that this contribution is subdominant 
at the scales of interest to us 

In the following calculations we mostly use the covariance 
matrix given in the middle row of Fig. 7.'^ 


5. Relation to the true spectrum 

Since masking in real space is equivalent to convolution in 
Fourier space, our measured power spectrum P is actually a 
convolution of the real spectrum P with a survey window (see 
e.g. FKP): 

E(k) = / ^20) 

where 

W(k) = J" d^rn(r)w(r) exp(/k ■ r) , (21) 

and the survey window |W(k)|^ is normalized as follows: 

( 22 ) 

The angle averaged survey window |W(k)p is plotted in Fig. 8. 
Here the core part of the window is well approximated by the 
functional form: 

\W{k)\^ =- \ -j, (a ^ 0.0030, b 0.0028), (23) 

i + (!) +(l) 

and asymptotic wings are close to the power law with spectral 
index -4. These approximations are shown with dashed lines in 

Since now we are also able to handle quasi-nonlinear mode-mode 
couplings. 

In calculating this contribution to the covariance matrix we have 
taken the best fit model parameters as obtained in Sec. 6. The small¬ 
ness of this term is caused by the high value of the parameter Mo i.e. 
majority of the “occupied” halos contain only one LRG. 

This matrix along with the power spectrum results in Fig. 5 is also 
given in a tabular form in Appendix G. 



wavenumber [hMpc'j 


Fig. 8 Isotropized survey window. Here the normalization is taken 
such that iiy(0)| = 1. Light gray stripe marks the region where the 
window is above 1% of its maximum value of 1. Dashed lines show 
approximations discussed in the text. 



Fig. 9 3D survey window embedded in a box with a side length of 
0.04 A Mpc^'. Here the isosurface corresponding to 1% of the maxi¬ 
mum value of the window is shown. Note the symmetry of the win¬ 
dow, |lV(k)p = |VT(-k)p, as expected when taking a modulus of the 
Fourier transform of a real 3D scalar function. 


Fig. 8. With the gray shaded stripe we have marked the scales 
where \W(k)\^ is above 1% of its maximum value. This stripe 
just serves as a rough guide to the effective width of the survey 
window and it is also shown in many of the following figures. 

Since the survey geometry of the analyzed SDSS LRG sam¬ 
ple is far from being spherically symmetric, an isotropized win¬ 
dow in Fig. 8 gives only a poor representation of the true 3D 
window, which is displayed as an isosurface corresponding to 
the isovalue of 0.01 in Fig. 9. 

In order to compare theoretical models to the measured 
power spectrum we have to take into account the smearing ef- 
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Fig. 10 Coupling kernels Al,(k) = K{ki, k) for the power spectrum bins 
kj shown in Fig. 5. Numerically evaluated kernels are shown with solid 
lines. The dashed lines correspond to the fitting functions given in 
Appendix D. 

fects caused by the survey window. Using Eq. (20) we can ex¬ 
press an isotropized power spectrum as: 

P{k)^ J'^p(k) = Jdk'k'^P(k')K(k',k), (24) 

where the coupling kernels 

K(k', k) = K(K ^ ^ ■ ^25) 

Numerically evaluated coupling kernels along with the analyt¬ 
ical approximations (see Appendix D) for the analyzed galaxy 
sample are presented in Fig. 10. Here the solid lines corre¬ 
spond to the numerical results and the dashed ones repre¬ 
sent an analytical approximation. We have used the notation 
Ki(k) = K(ki, k) where k, denote the central values of the power 
spectrum bins shown in Fig. 5. 

6. Model spectra 

It is well known that redshift space distortions and nonlinear ef¬ 
fects modify simple linear spectra. In order to treat these effects 
we make use of the very successful analytical model - the Halo 
Model. For a nice review we refer the reader to Cooray & Sheth 
(2002) (see also Seljak 2000). The details of the model we use 
are presented in Appendix C. The model introduces four free 
parameters: Miow, a, Mq and y. Here Miow is the lower cutoff 
of the halo mass i.e. below that mass halos are assumed to be 
“dark”, a and Mq are the parameters of the mean of the halo 
occupation distribution {N\M}, which gives the average num¬ 
ber of galaxies per halo with mass M. We take {N\M} to be a 
simple power law: 

( 26 ) 

We prefer to use “coupling kernels” instead of the more common 
“window functions” since the word “window” has already been used 
to mean the modulus square of the Fourier transform of the weighted 
survey volume. 


Fig. 11 Icr error contours for the free model parameters. Best fit pa¬ 
rameter values are marked with crosses. 


The last parameter, y, is the amplitude factor for the virial ve¬ 
locities of galaxies inside dark matter halos. One dimensional 
velocity dispersion of the galaxies inside a halo with mass M is 
taken to follow the scaling of the isothermal sphere model: 


cr — y 



(27) 


where /?vir is the virial radius of the halo. 

For the model fitting we have used Fevenberg-Marquardt 
method as described in Press et al. (1992) with modifications 
(described in Appendix E) that allow us to incorporate corre¬ 
lations between the data points. As the input data we take the 
power spectrum estimates given in Fig. 5. The covariance ma¬ 
trix used is the one shown in the middle row of Fig. 7. We also 
perform fits where we use one additional power spectrum bin 
on a larger scale (not shown in Fig. 5). All of this data is given 
in a tabular form in Appendix G. The transfer functions needed 
for the linear spectra are taken from Eisenstein & Hu (1998). 
There the authors also provide transfer function fits where the 
baryonic acoustic oscillations have been removed. We use these 
“smoothed out” transfer functions in order to assess the signifi¬ 
cance of the oscillatory features we see in the data. Throughout 
this paper we have kept cosmology fixed to the best fit WMAP 
“concordance” model (Spergel et al. 2003). The implications 
for the cosmology, and especially for the dark energy equation 
of state parameter, are planned to be worked out in the future 
paper. 

As the cosmology is kept fixed, we have only four free pa¬ 
rameters. In order to eliminate some of the degeneracies be¬ 
tween the parameters we have imposed one additional con¬ 
straint. Namely, we have demanded that the resulting number 
of galaxies should agree with the one that is observed with the 
relative error of 1% i.e. (51,763 + 518)'®. The resulting Icr 
error “ellipses” for the free parameters are shown in Fig. 11. 
The “ellipses” appear deformed since instead of Miow and Mq 

The 1(T Poisson error in this case would be 228. The large-scale 
structure amplifies the variability in the number of objects and a factor 
of a few increase above the Poissonian case seems to be reasonable. 
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we have fitted log(Miow) and log(Mo). With crosses we have 
marked the best fit values: Miow - 3 ■ Mq, a ^ 0.9, 

Mq ^ 1.4 ■ 10'‘*/i“'Mq and y ^ 0.7. The model spectra corre¬ 
sponding to these best fit parameters are shown in Figs. 4 and 
5. In both figures we have also given the simple linear spec¬ 
tra multiplied by the square of the bias parameter b - 1.95. 
In Fig. 4 we have additionally demonstrated the effect of the 
window convolution. There the dashed lines correspond to the 
unconvolved case. In Fig. 5 along with the “wiggly” spectra we 
have shown their “smoothed-out” counterparts. Using all the 16 
power spectrum bins (the 1st not shown in Fig. 5) plus an ad¬ 
ditional constraint on the total number of galaxies, resulting in 
17 - 4 = 13 independent degrees of freedom, we obtain val¬ 
ues of 8.8 and 19.9 for the “wiggly” and “smoothed”^® mod¬ 
els, respectively. So the models with oscillations are favored 
by 3.3cr over their “smoothed-out” counterparts.^' Since both 
models have the same number of free parameters, and if addi¬ 
tionally the assumption of Gaussianity is valid, the Bayesian 
approach should also give similar results. Actually, Bayesian 
results should favor “wiggly” models even more, since prior 
weight for these should probably be taken higher (assuming 
the knowledge of the other experimental results). 


7. Determination of the acoustic scale 

To measure the scale of the acoustic oscillations we divide the 
spectrum shown in Fig. 5 with the best fitting “smoothed” spec¬ 
trum. The result of this procedure is given in the upper panel of 
Fig. 12. There the solid line shows a cubic spline fit to the data 
points and the long-dashed line corresponds to the best fitting 
model spectrum also shown in Fig. 5. The above data is fitted 
with a parametric form: 



Fig. 12 Upper panel: Power spectrum from Fig. 5 divided by the best 
fitting “smoothed” spectrum. Solid line shows a cubic spline fit to 
the data points and long-dashed line corresponds to the best “wig¬ 
gly” model. The short-dashed line represents the most favorable fit 
from the parametric family of Eq. (28). Lower panel: Various input 
power spectra used to calculate the two-point correlation function. The 
dashed line is the cubic spline fit from the upper panel. The solid lines 
represent a transition sequence from the best fitting “wiggly” model 
to the best “smoothed” model. In each step we have erased more and 
more oscillatory features. For clarity slight vertical shifts have been 
introduced. 


f{x) = 1 -H Cl ■ sin(c 2 • x) exp 



(28) 


Again we use the Levenberg-Marquardt method with the 
data covariance matrix obtained from mock catalogs. After 
marginalizing over the other parameters we find the best fit¬ 
ting value of (105.4 + 23) Mpc for the parameter C 2 . 
The best fitting member of the parametric family in Eq. (28) 
is shown with short-dashed lines in the upper panel of Fig. 12. 
Using FKP covariance matrix instead gives an acoustic scale of 
(105.4 +2.8)Mpc. 


The best fit a, Mq and y for the “smoothed” models differ 
slightly from the values quoted above for the “wiggly” spectra. 

Dropping the first power spectrum bin the obtained values are 
5.0 and 16.5. 5.0 is an anomalously low value of for 12 degrees 
of freedom. (One would expect )^ - 12 ± 5.) In fact, if we would 
have used the simple FKP covariance matrix instead of the one ob¬ 
tained from the mock catalogs, the resulting values would be even 
lower: 2.9 and 8.5, respectively. This might hint that the 2nd order 
Lagrangian approach, although very successful, might still have prob¬ 
lems of capturing some extra mode-mode couplings. 

Here and in the following all the errors refer to the l-cr level. 
Values for the other parameters are as follows: ci = (4.9 ± 2.1) • 10“^, 
C3 = (0.176 ±0.023) AMpc-',C4 = (7 ± >’). 


The sinusoidal modulation in the power spectrum is a pure 
consequence of the adiabatic initial conditions. By relaxing this 
assumption and fitting with a more general functional form: 


f(x) = 1 + Cl • sin(c 2 • X + C 3 ) exp 



(29) 


instead, we get the following value for the acoustic scale: 
(103.0 + 7.6)/i“' Mpc. In case of the FKP covariance matrix 
the corresponding value is (103.1 + 9.1) Mpc. 

Eisenstein et al. (2005) determine various distance scales 
(like which is a certain mixture of the comoving distances 
along and perpendicular to the line of sight (see their Eq.(2))) 
and their ratios, using SDSS LRGs in combination with the 
constraints from other cosmological sources. The typical rel¬ 
ative accuracy of these measurements is ~ 4%, which might 
seem to be significantly poorer than the accuracy of the acous¬ 
tic scale measurement, (105.4 + 2.3) Mpc i.e. ~ 2%, pre¬ 
sented in this paper. This apparent inconsistency can be at¬ 
tributed to the fact that in our analysis, as stated above, we 
have kept the cosmology fixed to the WMAP “concordance” 
model, whereas Eisenstein et al. (2005) estimates include the 
extra uncertainties due to the imperfect knowledge of the vari¬ 
ous cosmological parameters. Of course, the given length of the 
acoustic scale, (105.4 + 2.3) Mpc, can be easily transformed 
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in order to accommodate other preferences for the background 
cosmology. We also note that the use of the parametric form in 
Eq. (28) might be too restrictive, since the acoustic modulation 
in the case of adiabatic models can be only approximately de¬ 
scribed as a damped sinusoidal wave (Eisenstein & Hu 1998). 
Eor this reason the given sound horizon constraint should not 
be used in cosmological parameter studies. Instead one should 
directly use the measured power spectrum in combination with 
the parametrized models that are physically well motivated. 


8. Correlation function analysis 


We determine the two-point correlation function of the SDSS 
LRGs using the edge-corrected estimator given by Landy & 
Szalay (1993): 


^(0 = 


DD - 2DR + RR 


(30) 


which has minimal variance for a Poisson process. Here DD, 
DR and RR represent the respective normalized data-data, data- 
random and random-random pair counts in a given distance 
range. Random catalogs were generated with 25 times the num¬ 
ber of objects in the main catalog. We calculated correlation 
function for 10 bins (r,, i - 1... 10) in the pair distance range 
of 60... 160 Mpc. The errors were estimated by a “jack¬ 
knife” technique. Eor this purpose we divided the full sample 
into 22 separate regions each containing ~ 2350 galaxies (see 
Eig. 3). The two-point function was calculated 22 times, each 
time omitting one of the regions. Denoting the resulting esti¬ 
mates as (y = 1... 22), the “jackknife” estimate for the 

variance reads as (see e.g. Lupton 1993): 

^31) 

y=l 

1 ^ 

^32) 

J=l 

where in our case N - 22. The results of this calculation are 
presented in the left panel of Eig. 13. With the crosses and 
dashed-line errorbars we have also shown the two-point func¬ 
tion as determined by Eisenstein et al. (2005). We see that in 
general our results agree reasonably well with their calcula¬ 
tions. 

It would be interesting to study how the oscillations in the 
observed power spectrum transform into the peak in the two- 
point correlation function seen at the scale of ~ llO/i ' Mpc. 
Eor this purpose we use the cubic spline fit shown in Eig. 12 
and extend it outside of the observed range by smoothly join¬ 
ing it to the power spectrum of the best fitting “smoothed-out” 
model. The correlation function is now simply calculated as the 
Eourier transform of the power spectrum. The resulting cor¬ 
relation function is plotted with a dashed line in the right panel 

To be precise, in redshift space the two-point correlation func¬ 
tion and power spectrum are not anymore exact Fourier transforms of 
each other. Nevertheless, we think that this simplified exercise is still 
useful. Also, as the correlation function estimator in Eq. (30) is an 
edge-corrected estimator, we use an unconvolved model spectra here. 



r[h''Mpc] r[h‘^Mpc] 

Fig. 13 Left panel: Two-point correlation functions as determined in 
this paper (circles with solid lines) and hy Eisenstein et al. (2005). 
Right panel: Correlation functions corresponding to the models shown 
in the lower panel of Eig. 12 in comparison to the one obtained directly 
from the data. Here all the data points have been lowered by 0.0035. 


of Fig. 13. To study the significance of the oscillatory features 
in the power spectrum in relation to the observed peak in the 
correlation function, we have calculated correlation functions 
for several models that have oscillations “switched off” at var¬ 
ious scales. The spectra of these models are shown with solid 
lines in the lower panel of Fig. 12, where for the sake of clarity 
we have introduced slight vertical shifts between the curves, so 
that the scales where the transition to the featureless spectrum 
takes place, are easily visible. The corresponding correlation 
functions are given with solid lines in the right hand panel of 
Fig. 13. As expected, we see how the peak in the correlation 
function is getting broader and also decreasing in amplitude as 
we erase more and more features in the power spectrum. This 
clearly demonstrates the importance of many of the up-downs 
in the power spectrum to produce a relatively sharp feature in 
the two-point correlation function. 

In order to achieve good agreement we have lowered all the 
data points by 0.0035 in the right hand panel of Fig. 13. Similar 
shifts were also suggested in Eisenstein et al. (2005) in order 
to get better match to the theoretical models. A 0.0035 shift in 
f translates to the 0.175% shift in the mean density. Thus, if 
one wishes to determine the amplitude of the correlation func¬ 
tion correctly at those large scales, one has to determine the 
survey selection function with a very high precision, which in 
practice is very difficult to achieve. By using model spectra that 
have more large scale power than the “concordance” cosmol¬ 
ogy predicts (as might be suggestive from Fig. 4), we are in fact 
able to match the amplitude of the correlation function without 
any additional vertical shifts. Here we try to avoid making any 
definite conclusions. The behavior of the power spectrum on 
the largest scales is an extremely interesting topic on its own 
and there exist much better methods than the direct Fourier ap¬ 
proach to investigate these issues (see e.g. Tegmark et al. 1998). 
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Fig. 14 The comparison of spectra from different surveys. 



k [hMpc'h 

Fig. 15 The same as Fig. 14 with the errorbars omitted. 

9. Comparison with the other surveys 

In this section we compare our power spectrum measurements 
with the ones obtained by Percival et al. (2001) and Cole et al. 
(2005) for the 2dF redshift survey and by Tegmark et al. (2004) 
for the SDSS main galaxy sample. The results of this compari¬ 
son are provided in Figs. 14 and 15. Since Fig. 14 is extremely 
busy, we have also given a variant of it where we have omitted 
the errorbars. The amplitudes of the SDSS main and 2dF spec¬ 
tra have been freely adjusted to match the clustering strength 
of the SDSS LRGs. The corresponding bias parameters with 
respect to the SDSS LRGs are 0.53, 0.61 and 0.50 for the 2dF 
sample analyzed by Percival et al. (2001), for the one analyzed 
by Cole et al. (2005), and for the SDSS main sample, respec¬ 
tively. Percival et al. (2001) also provide power spectrum mea¬ 
surements for k > 0.15 hMpc * but without errorbars. These 
small-scale measurements along with our SDSS LRG results 
are shown with solid lines in Fig. 14. 

In general the shapes of the spectra agree remarkably well. 
Of course one has to keep in mind that here, with the only ex¬ 
ception of Tegmark et al. (2004) results, the power spectrum 


bins are highly correlated. Also Tegmark et al. (2004) mea¬ 
surements are corrected for the redshift space distortions. 

10. Discussion and Conciusions 

In this paper we have calculated the redshift-space power spec¬ 
trum of the SDSS DR4 LRG sample, hnding evidence for a 
series of acoustic features down to the scales of ~ 0.2 h Mpc“*. 
It turns out that models with the baryonic oscillations are fa¬ 
vored by 3.3cr over their “smoothed-out” counterparts with¬ 
out any oscillatory behavior. Using the obtained power spec¬ 
trum we predict the shape of the spatial two-point correlation 
function, which agrees very well with the one obtained directly 
from the data. Also, the directly calculated correlation func¬ 
tion is consistent with the results obtained by Eisenstein et al. 
(2005). We have made no attempts to put constraints on the cos¬ 
mological parameters, rather we have assumed in our analysis 
the “concordance” cosmological model. The derived acoustic 
scale (105.4+2.3) Mpc agrees well with the best-ht WMAP 
“concordance” model prediction of ^ 106.5 /j * Mpc. 

The existence of the baryonic features in the galaxy power 
spectrum is very important, allowing one (in principle) to ob¬ 
tain Hubble parameter H and angular diameter distance as 
a function of redshift, this way opening up a possibility to con¬ 
strain properties of the dark energy (Hu & Haiman 2003). The 
currently existing biggest redshift surveys, which are still quite 
shallow, do not yet provide enough information to carry out 
this project fully. On the other hand, it is extremely encour¬ 
aging that even with the current generation of redshift surveys 
we are already able to see the traces of acoustic oscillations 
in the galaxy power spectrum, showing the great promise for 
the dedicated future surveys like K.A.O.S. We have seen that 
acoustic features seem to survive at mildly nonlinear scales 
{k > 0.1 /iMpc^'), which is in agreement with the results of 
the recent N-body simulations (Springel et al. 2005; Seo & 
Eisenstein 2005). In order to fully exploit available information 
one needs a complete understanding of how nonlinear effects 
influence these features. Nonlinear bias and redshift space dis¬ 
tortions also add extra complications. In general redshift-space 
distortions, biasing and nonlinear evolution do not create any 
oscillatory modulation in the power spectrum and so acoustic 
features should be readily observable. So far there have been 
only a few works studying these important issues (e.g. Springel 
et al. (2005); Seo & Eisenstein (2005); White (2005)) and prob¬ 
ably it is fair to say that currently we really do not have a full 
theoretical description of them. In our paper we have mod¬ 
eled the above mentioned effects using the results from the 2nd 
order Lagrangian perturbation theory in combination with the 
Halo Model. Although these models are very successful in cap¬ 
turing many important aspects of the structure formation, one 
has to keep in mind that they are still approximations. 

The bare existence of the baryonic oscillations in the galaxy 
power spectrum tells us something important about the under¬ 
lying cosmological model and the mechanism of the structure 
formation. Eirst, it confirms the generic picture of the gravi¬ 
tational instability theory where the structure in the Universe 
is believed to be formed by the gravitational amplihcation 
of the small perturbations layed down in the early Universe. 
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Under the linear gravitational evolution all the density fluctu¬ 
ation modes evolve independently i.e. all the features in the 
power spectrum will be preserved. And certainly, we are able 
to identify features in the low redshift galaxy power spectrum 
that correspond to the fluctuations seen in the CMB angular 
power spectrum (which probes redshifts z ~ 1100), providing 
strong support for the above described standard picture of the 
structure formation. Actually, we can also probe scales that are 
unaccessible for the CMB studies due to the strong damping ef¬ 
fects and steeply rising influence of the secondary anisotropies, 
reaching effectively the wavenumbers that correspond to the 
6th-7th peak in the CMB angular power spectrum. Second, the 
ability to observe baryonic features in the low redshift galaxy 
power spectrum demands rather high baryonic to total matter 
density ratio. In Blanchard et al. (2003) it has been shown that 
it is possible to fit a large body of observational data with an 
Einstein-de Sitter type model if one adopts low value for the 
Hubble parameter and relaxes the usual assumptions about the 
single power law initial spectrum. In the light of the results ob¬ 
tained in our paper these models are certainly disfavored due 
to the fact that the high dark matter density completely damps 
the baryonic features. And finally, purely baryonic models are 
also ruled out since for them the expected acoustic scale would 
be roughly two times larger than observed here . So the data 
seems to demand a weakly interacting nonrelativistic matter 
component and all the models that try to replace this dark mat¬ 
ter component with something else e.g. modifying the laws of 
gravity might have severe difficulties to fit these new observa¬ 
tional constraints. 
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Fig.A.l Power spectra of galaxy clusters from the z = 0 Hubble 
Volume simulation box. The solid line represents the spectrum for the 
full sample of 1,560,995 clusters. The circles with errorbars denote 
the recovered spectrum from the 50 mock catalogs having similar se¬ 
lection effects to the analyzed SDSS LRG sample. The dotted lines 
demonstrate the convolution effect of the survey window on the best 
fitting model spectrum. 
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Appendix A: Test problem 

Here we present one test for our power spectrum calculation soft¬ 
ware As the input we use the z = 0 cluster catalog of the VIRGO 
Hubble Volume simulations^®, which covers the comoving volume of 
3000Mpc^ and contains 1,560,995 clusters above the mass limit 
of 6.75 • 10 '^/!“'Mq. The average bias parameter of this catalog is 
b - 1.9, which is comparable to the SDSS LRG value of = 1.95. The 
power spectrum of the full sample is shown in Fig. A.l with a solid 
line. Here for clarity we have not shown the errorbars, which are rather 
small for a sample of that size. Out of the full sample we generate 50 
mock catalogs that have the same radial and angular selection func¬ 
tions as the SDSS LRG sample analyzed in this paper (see Figs. 2 and 
3). The mean number of objects in the resulting catalogs is ~ 18,500 

i.e. the number density is roughly one third of the spatial density of the 
SDSS LRGs. Observer’s location and pointing angles are taken ran¬ 
domly for each of the catalogs. The mean recovered power spectrum 
with Icr errorbars is shown in Fig. A.l. We see that the power spectrum 
of the underlying sample is recovered very well. On the largest scales 
there are some deviations, which can be explained as being caused by 
the smearing effect of the survey window. This is demonstrated by the 
dotted lines, where the lower/upper curve corresponds to the model 
spectrum with/without survey window convolution applied. 

Further tests can be found in Hiitsi (2005) 

http://www.mpa-garching.mpg.de/Virgo/ 


Appendix B: Mock catalogs 

We build mock catalogs for the SDSS LRG by a 3 step procedure: 

1. Generation of the density field using an optimized 2nd order 
Lagrangian perturbation calculation (2LPT). 

2. Poisson sampling of the generated density field with the intensity 
of the process adjusted so, as to end up with a galaxy sample that 
has a clustering strength enhanced by a factor with respect to 
the underlying field, and a number density equal to the observed 
LRG sample density at the minimal used redshift of 0.16 (see Fig. 
2 ). 

3. Extraction of the final catalog by applying the radial and angular 
selection function as given in Figs. 2 and 3, respectively. 

In contrast to the Eulerian perturbation theory, where one does a 
perturbative expansion of the density contrast field, Lagrangian ap¬ 
proach considers an expansion of the particle trajectories (see e.g. 
Buchert & Ehlers 1993; Bouchet et al. 1995; Sahni & Coles 1995; 
Bernardeau et al. 2002). Here the central quantity is the displace¬ 
ment field 'l'(q), which relates particle’s initial comoving position 
(Lagrangian position) q to its final Eulerian position x: 

x = q-H'P(q). (B.l) 

It turns out that due to the decay of the rotational perturbation modes in 
the expanding Universe each order of the perturbation theory displace¬ 
ment field separates into a time-dependent and a Lagrangian position 
dependent factors (Ehlers & Buchert 1997). The position dependent 
part, due to its irrotational nature can be given as a gradient of a scalar 
potential. As a result, one can expand the displacement field as fol¬ 
lows: 

Tiq) = -t . (B.2) 

Here the 1st term describes the classical Zeldovich approximation 
(Zel’Dovich 1970). The time independent potentials 0 ^) ^( 2 ) ^j-g 

found from the Poisson equations: 

A 0 <‘>(q) = -d(q) (B.3) 

and 

A 0 “>(q) = 5 E Z ’ (B.4) 

' / 

where denotes the partial derivative with respect to the Lagrangian 
coordinate d(q) is the initial density contrast. We generate d(q) us¬ 
ing the standard Zeldovich approximation on a regular cubical grid. 

Di in Eq. (B.2) is the linear growth factor. The second-order 
growth factor D 2 for flat models with a cosmological constant is to 
a good precision approximated as (Bouchet et al. 1995): 

Dj ^ ^B.5) 

According to Eq. (B.l) and (B.2) the peculiar velocity field is given 
as: 

v = -H . (B, 6 ) 

Here H = ^ and . For flat models with a cosmological 

constant logarithmic derivatives of the growth factors can be approxi¬ 
mated as (Bouchet et al. 1995): 

/i - /2 - 2D®/‘i . (B.7) 

Lagrangian perturbative approach works fine up to the 1st shell¬ 
crossing. After that the formed caustic structures will start to be wiped 
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Fig. B.l A 25/!“' Mpc thick slice through a 1280/!“' Mpc computational box. A gray scale image represents the underlying density field 
obtained by the optimized 2LPT approach. White dots mark the positions of the “galaxies” generated by the Poisson sampler. 


out, since the particles just keep on moving without noticing the grav¬ 
itational pull of the dense sheets/filaments. It is possible to cure this 
problem significantly by filtering out the small-scale Fourier modes. 
This is what is meant by the “optimization”. The method applied to the 
1st order Lagrangian perturbation calculation is known as the trun¬ 
cated Zeldovich approximation (e.g. Coles et al. 1993; Melott et al. 
1994; Weiss et al. 1996). Weiss et al. (1996) suggest to remove the 
small-scale power by applying a Gaussian k-space filter with a charac¬ 
teristic smoothing scale kj., to the initial density field. Thus the power 
spectrum of the filtered field is given by: 


Poptimized(k) — P(k) exp 



(B.8) 


They recommend the value kgs - 1.2k„;, where the nonlinearity scale 
k„i is defined as: 


(2ny 


^nl 

f 


d^kP(k) = 1. 


0 


(B.9) 


Although they studied only models with fl„ = 1, it has been later 
shown by Hamana (1998) that this “recipe” performs well for arbitrary 
Friedmann-Lemaitre-Robertson-Walker models. 

In our calculations we assume the WMAP “concordance” cos¬ 
mology (Spergel et al. 2003). Linear power spectrum is taken from 
Eisenstein & Hu (1998). We build 2LPT density field on a 256^-grid 
with 5 /!“' Mpc cell size using the same number of particles as the 
number of grid cells. Four copies of this box are combined to form 
a bigger 2560 x 2560 x 1280/i“^ Mpc^ volume. Out of that big box a 
sample of “galaxies” is selected with a radial number density as given 
in Fig. 2 and with an angular mask presented in Fig. 3. The parameters 
of the Poisson sampler are tuned to give a sample with a bias param¬ 
eter h 2 in agreement with the observed value for the SDSS LRG 

Due to the rather big cell size the truncation of the initial spectrum 
has a rather mild effect. 

We use a simple model where the intensity of the inhomogeneous 
Poisson process is linearly related to the underlying density field. 
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and the 2-halo term: 


P^\k) = I r" + + \rn Pijk). 


(C.4) 


Here: 


•Rja = ka- ^ = 


-0r erf(a) 


(C.5) 




(N\M) 

dM n(M)b(M)^-^‘Ri(kcr)ug(k\M), (C.6) 


/ 


Ty = f- dMn(M)b(M)'Ri(ko-)u(k\M). 


(C.7) 


In the above expressions n(M) is the mass function and b(M) halo bias 
parameter. We calculate them using a prescription by Sheth & Tormen 
(1999) and Sheth et al. (2001). h represents the mean number density 
of galaxies: 


Fig. B.2 The power spectrum of ~ 350,000 “galaxies” from the simu¬ 
lation box shown in Fig. B. 1. The solid line shows the linearly evolved 
input spectrum multiplied by the square of the bias parameter b = 2.0. 


sample. The redshift-space catalog is built by altering the radial dis¬ 
tances of the “galaxies” by \r/Ho, where Vr is the radial component of 
the peculiar velocity field (see Eq. (B.6)) and Hq = 100 Akm/s/Mpc. 

In Fig. B.l we show a 25 Mpc thick slice through a box with 
1280Mpc side length. The underlying density field is presented 
as a gray scale image with white dots marking the positions of the 
“galaxies”. The power spectrum of the sample of ~ 350,000 “galax¬ 
ies” is shown in Fig. B.2 We see that the shape of the spectrum is 
in good agreement with the linearly evolved power spectrum up to the 
scales of k ~ 0.5 iiMpc^'. 

This approach gives us a “galaxy” sample that has rather realis¬ 
tic large-scale clustering properties. In terms of the Halo Model (see 
Appendix C) one can say that halo-halo clustering term is properly ac¬ 
counted for. 2LPT also gives reasonably accurate higher order correla¬ 
tions on quasi-nonlinear scales (e.g. Bouchet et al. 1995; Scoccimarro 
& Sheth 2002). 

Appendix C: Power spectrum from the halo model 


n = 


f 


dMn(M)(NIM}. 


(C.8) 


We take the mean of the halo occupation distribution in the following 
form: 


(N\M) = 



(C.9) 


where Mg and a are free parameters. The second moment is chosen as 
(see Cooray 2004): 


(N(N - 1)|M) 

m) 




I l‘°g( lonr-iMe ) 


(C.IO) 


if M < 10*3 
otherwise. 


(C.ll) 


/ in Eq. (C.7) denotes the logarithmic derivative of the linear growth 
factor: / = . u{k\M) and Ug{k\M) are the normalized Eourier 

transforms of the dark matter and galaxy density distributions within a 
halo of mass M. In our calculations we take both of these distributions 
given by the NEW profile (Navarro et al. 1997) and the concentration 
parameter c{M) is taken from Bullock et al. (2001). One dimensional 
velocity dispersion of the galaxies inside a halo with mass M is taken 
to follow the scaling of the isothermal sphere model: 


The halo model description of the spatial clustering of galaxies is a 
development of the original idea by Neyman & Scott (1952), where 
one describes the correlations of the total point set as arising from the 
two separate terms: (i) 1-halo term, that describes the correlations of 
galaxies populating the same halo, (ii) 2-halo term, which takes into 
account correlations of the galaxies occupying different halos. For a 
thorough review see Cooray & Sheth (2002). Here we briefly give the 
results that are relevant to the current paper (see Seljak 2001; Cooray 
2004). 

The power spectrum of galaxies in redshift space can be given as: 
P{k) = P^’'(k) + (C.l) 

where the 1-halo term: 

= J dM n(M) ~ 1ip(kcr)\Ug(k\M)\'’ , (C.2) 

r 1 if {N(N - DIM) <1 

^ \ 2 if {N(N - DIM) > 1 3'— 

Here as in the previous figure the SDSS LRG selection functions 
are not applied yet. 


cr - y 



(C.l 2) 


where Ryir is the virial radius of the halo and y is a free parameter. 

After specifying the background cosmology the above described 
model has four free parameters: Mg, a (Eq. (C.9)), cr (Eq. (C.12)) and 
Mio„. The last parameter Mio„ represents the lower boundary of the 
mass integration i.e. halos with masses below Mio„ are assumed to be 
“dark”. 

One can also use the halo model to estimate nonlinear contribu¬ 
tions to the power spectrum covariance matrix. The additional term to 
the covariance matrix (i, y-denote power spectrum bins) arising 
from the parallelogram configurations of the trispectrum is given by 
(Cooray 2004): 


Tij 1 r d?k Pd^k r 
= — = - - - dMn(M)- 

V Vj V, J Vj J 


r'NL _ 


{N(N - 1)(N - 2)(N ■ 


J 

■ 3)|M) 


Ug{k\Mf\Ug{kj\Mf 


(C.13) 


3° Here only the contribution due to the 1-halo term is given. 
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f denotes an integral over a fe-space shell centered at wavenumber k/ 

i 

with a volume Vi = AnkjAk. The 4th moment of the halo occupation 
distribution is taken as: 

{N{N - 1)(N - 2)(N - 3)|M) = - l] • 

[3^^(M) - 2 ] <iV|M>^ (C.14) 

where j8(M) and (N\M) are given in Eqs. (C.ll) and (C.9) above. 
Performing calculations in redshift space a factor of 'Rp(kcr) (see Eq. 
(C.5)) must also be included in Eq. (C.13). 


Appendix D: Fitting formuiae for the coupiing 
kerneis 


In this appendix we provide analytical fitting formulae for the cou¬ 
pling kernels K(k, k') in Eq. (24) . The analytic form is motivated by 

the fact that the angle averaged survey window |VT(k)P (see Fig. 8) can 
be reasonably well approximated by the analytical form: 


\w{kf = m = 


1 



(D.l) 


Now assuming that |iy(k)p is isotropic (which certainly is not the case 
as seen from Fig. 9), we can find the coupling kernels K(k, k') as: 


K(k,k') = C- 





k+k' 

9,71^ r 

C-— J /(x)xdx.(D.2) 

lk-k'\ 


For f(k) given by Eq. (D. 1) the integral in Eq. (D.2) and the normaliza¬ 
tion constant C can be found analytically. The kernels are normalized 
such that 


f 


K(k,k')k'^dk' = 1 


(D.3) 


is satisfied. 

Depending on the values of a and b there are two different solu¬ 
tions. 


1. b* > 40*: 

g(k - k') 
g(k + k') 

where 

/! + b‘^ + 2a^x^ 

" II-b*- 2a2x2 ’ 

11 = b^ VM - 4a‘^ 

and the normalization constant: 

C =_ ^ . 

+ 11 - x/b‘* - iij 

2. b* < 4a*-. 

K(k, k' ) = ^ [g{k -t k') - g(k - k' )] , 


K{k, k') = In 
kk' 


where 
g(x) = arctan I 


b* + 2a^x^ 


W^4a*-b*l 
and the normalization constant: 
1 

C = 




(D.4) 


(D.5) 

(D.6) 

(D.7) 


(D.8) 


(D.9) 


(D.IO) 


To avoid confusion we do not call them window functions since 
the word “window” has been already used to mean the Fourier trans¬ 
form of the survey volume. 


Although the isotropy assumption is certainly not correct, the above 
parametric family provides a very good fit to the numerically evaluated 
kernels as seen in Fig. 10. The best fitting a and b for the analyzed 
SDSS LRG sample are 0.00457 and 0.00475, respectively. 


Appendix E: Nonlinear model fitting. Correlated 
data 


We find the best fitting parameters for the nonlinear model by mini¬ 
mizing which in the case of Gaussian errors is equivalent to finding 
the maximum likelihood solution. For this purpose we use Levenberg- 
Marquardt method as described in Press et al. (1992), where it was 
assumed that data values are uncorrelated. Since we are interested in 
the case with correlated errors, we have to make slight modifications 
to their implementation of the algorithm. 

Using their notation, is now calculated as: 

N N 

X^(a) = 2 ■ ‘“d' ■ 

,=1 ;=1 


and the quantities and as follows: 




,=i j=i 


j dy(Xj-,a) 
' dai. 


(E.2) 


au 


dyjxj-, a) [ d-yjxf a) 
dai ■ '> ■ da, 

1=1 j=i 


(E.3) 


In the above relations Cij represents the data covariance matrix. 


Appendix F: Goodness of fit. Correlated Gaussian 
data 

Under the assumption that statistical fluctuations Ay, = y,- - y(x,-; a) 
(i = 1... A) in Eq. (E.l) are Gaussian distributed, with covariance 
matrix C,y, one can easily derive probability density function (pdf) 
for the quantity and thus open up a way to estimate the goodness 
of fit. goodness-of-fit estimator is usually exploited in the case of 
independent Gaussian variables. Here we show that calculating^^ for 
the correlated Gaussian data as given in Eq. (E.l), one obtains the 
same result that is well known for the independently distributed case. 
According to our assumption Ay is Gaussian distributed: 

/Ay(Ay) = , ^ - exp (-^Ay^ • C^' • Ay]. (E.l) 

V(27r)'VdetC \ 2 j 

The conditional pdf ofx^ given Ay: 

/UlAyOr'lAy) = s{x^- Ay^ • C-‘ • Ay), (F.2) 

and so the pdf for;y^ can be written as: 

f, 2 (x^) - , ^ f d'^Ay expf-^Ay^ • C-' • Ay] • 

V(2;r)~detC J \ 2 j 

6{x^-Ay^ -C-' -Ay). (F.3) 

Now we define a new set of variables: 

Ay' = • Ay, (F.4) 

where L is the lower triangular matrix appearing in the Cholesky de¬ 
composition of C“‘: 

C“‘=L-L^. (F.5) 



The metric in the new frame is an identity matrix. 



bin # 

ki [h Mpc ^ ] 


A/>[A-5Mpc3] 

1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

1 

0.148e-l 

0.133e6 

0.228e5 

0.518e9 

0.358e8 

0.762e7 

0.521e7 

0.337e7 

0.161e7 

0.189e7 

0.987e6 

0.128e7 

0.106e7 

0.400e6 

0.246e6 

0.204e6 

0.114e6 

0.106e6 

0.184e6 

2 

0.315e-l 

0.712e5 

0.461e4 

0.358e8 

0.212e8 

0.267e7 

0.122e7 

0.433e6 

0.402e6 

0.260e6 

0.184e6 

0.198e6 

0.165e6 

0.146e6 

0.124e6 

0.958e5 

0.854e5 

0.420e5 

0.317e5 

3 

0.504e-l 

0.500e5 

0.224e4 

0.762e7 

0.267e7 

0.500e7 

0.103e7 

0.362e6 

0.316e6 

0.192e6 

0.170e6 

0.810e5 

0.702e5 

0.515e5 

0.467e5 

0.389e5 

0.285e5 

0.138e5 

0.151e5 

4 

0.697e-l 

0.394e5 

0.140e4 

0.521e7 

0.122e7 

0.103e7 

0.196e7 

0.399e6 

0.204e6 

0.159e6 

0.142e6 

0.117e6 

0.986e5 

0.839e5 

0.614e5 

0.652e5 

0.492e5 

0.152e5 

0.215e5 

5 

0.891e-l 

0.286e5 

0.816e3 

0.337e7 

0.433e6 

0.362e6 

0.399e6 

0.667e6 

0.160e6 

0.112e6 

0.110e6 

0.788e5 

0.654e5 

0.537e5 

0.647e5 

0.564e5 

0.477e5 

0.266e5 

0.298e5 

6 

0.109e0 

0.213e5 

0.550e3 

0.161e7 

0.402e6 

0.316e6 

0.204e6 

0.160e6 

0.303e6 

0.113e6 

0.798e5 

0.614e5 

0.520e5 

0.416e5 

0.427e5 

0.422e5 

0.325e5 

0.233e5 

0.265e5 

7 

0.128e0 

0.185e5 

0.462e3 

0.189e7 

0.260e6 

0.192e6 

0.159e6 

0.112e6 

0.113e6 

0.213e6 

0.817e5 

0.594e5 

0.436e5 

0.384e5 

0.284e5 

0.299e5 

0.256e5 

0.215e5 

0.217e5 

8 

0.148e0 

0.147e5 

0.364e3 

0.987e6 

0.184e6 

0.170e6 

0.142e6 

0.110e6 

0.798e5 

0.817e5 

0.133e6 

0.533e5 

0.366e5 

0.312e5 

0.269e5 

0.217e5 

0.188e5 

0.181e5 

0.129e5 

9 

0.167e0 

0.119e5 

0.295e3 

0.128e7 

0.198e6 

0.810e5 

0.117e6 

0.788e5 

0.614e5 

0.594e5 

0.533e5 

0.871e5 

0.384e5 

0.266e5 

0.242e5 

0.217e5 

0.188e5 

0.128e5 

0.143e5 

10 

0.187e0 

0.103e5 

0.251e3 

0.106e7 

0.165e6 

0.702e5 

0.986e5 

0.654e5 

0.520e5 

0.436e5 

0.366e5 

0.384e5 

0.631e5 

0.254e5 

0.200e5 

0.169e5 

0.136e5 

0.927e4 

0.100e5 

11 

0.206e0 

0.869e4 

0.224e3 

0.400e6 

0.146e6 

0.515e5 

0.839e5 

0.537e5 

0.416e5 

0.384e5 

0.312e5 

0.266e5 

0.254e5 

0.504e5 

0.237e5 

0.181e5 

0.158e5 

0.114e5 

0.123e5 

12 

0.226e0 

0.759e4 

0.210e3 

0.246e6 

0.124e6 

0.467e5 

0.614e5 

0.647e5 

0.427e5 

0.284e5 

0.269e5 

0.242e5 

0.200e5 

0.237e5 

0.44 le5 

0.217e5 

0.174e5 

0.148e5 

0.141e5 

13 

0.246e0 

0.654e4 

0.199e3 

0.204e6 

0.958e5 

0.389e5 

0.652e5 

0.564e5 

0.422e5 

0.299e5 

0.217e5 

0.217e5 

0.169e5 

0.181e5 

0.217e5 

0.396e5 

0.191e5 

0.133e5 

0.146e5 

14 

0.265e0 

0.599e4 

0.182e3 

0.114e6 

0.854e5 

0.285e5 

0.492e5 

0.477e5 

0.325e5 

0.256e5 

0.188e5 

0.188e5 

0.136e5 

0.158e5 

0.174e5 

0.191e5 

0.330e5 

0.152e5 

0.120e5 

15 

0.285e0 

0.530e4 

0.162e3 

0.106e6 

0.420e5 

0.138e5 

0.152e5 

0.266e5 

0.233e5 

0.215e5 

0.181e5 

0.128e5 

0.927e4 

0.114e5 

0.148e5 

0.133e5 

0.152e5 

0.261e5 

0.155e5 

16 

0.305e0 

0.479e4 

0.164e3 

0.184e6 

0.317e5 

0.151e5 

0.215e5 

0.298e5 

0.265e5 

0.217e5 

0.129e5 

0.143e5 

O.lOOeS 

0.123e5 

0.141e5 

0.146e5 

0.120e5 

0.155e5 

0.271e5 
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